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Abstract 

GZKFast is a low-cost astrophysical event generator designed to simulate photohadron processes 
resulting from ultra high energy cosmic ray fluxes. GZKFast is an easy to use event generator which 
specifically addresses issues relevant to cosmogenic neutrinos and the ultra high energy (UHE) neutrino 
spectrum. GZKFast injects UHE particles into a simulated cosmic microwave background (CMB) using 
a Monte Carlo approach. The interaction of each particle is simulated and the resulting detectable 
events are spooled to log files. Although GZKFast uses a simplified particle physics model and limited 
accelerator data it generates event distributions which are comparable with more thorough simulations 
such as the event generator SOPHIA [T|. 



1 Advertisement 

GZKFast is a capable event generator which pro- 
duces results comparable with applications that in- 
corporate far more substantial particle physics kine- 
matics and accelerator data. 

GZKFast provides... 

• A full featured simulation with rich run time 
configuration. 

• Monte Carlo event generation. 

• Management of individual astrophysical 
sources. 

• Individual event simulation and tracking in- 
cluding relevant particle kinematics. 

• Modular C++ design. 

• Histogram and extensive tabular output. 

• Practical standalone approach designed for 
ease of use, portability, and redistribution. 

• Reusable program library for easy incorpora- 
tion into new projects. 

• GNU Library General Public License. 



2 Introduction 

This year the "GZK Effect" of Greisen [2J and Zat- 
sepin and Kuz'min [3] is celebrating its fortieth 
anniversary. These authors independently realized 
that a uniform cosmic microwave background would 
provide an optically thick background for the high- 
est energy cosmic rays. Cosmic ray protons of suffi- 
ciently large energy would be likely to interact with 
CMB photons over suitably large distances. 

Accelerator experiments predict that the 
important channels for such interactions would be 
direct photopion production and the resonance chan- 
nels for the A g]. 

p + 7 — > n + 7r + 
p + 7^A + ^n + 7r + (1) 

. . . -> p + 7T° 

The stable daughters of these processes may 
be observable in Earth based detectors. Therefore 
we wish to predict the expected fluxes of these par- 
ticles. 

Four momentum invariance of these reac- 
tions implies that in the proton rest frame [5]: 



s = Mt 



where M p is the proton mass and E 7 is the 
energy of the photon in the lab frame. Therefore 
the measured cross section for photo-hadron inter- 
actions, er(s), determines the likelihood for proton 
attenuation. The expected attenuation length can 
be characterized by a mean free path, A. 
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The mean free path also depends on the av- 
erage density of CMB photons n avg . The average 
density may be found by assuming the CMB is char- 
acterized by a black body spectrum and integrating 
the Bose distribution. Including redshift, the aver- 
age density of the CMB is given by @ . 
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Here is the zeta function defined in sta- 
tistical mechanics texts and T is the temperature of 
the microwave background today. Using accepted 
values, at present this quantity is: 

n aV g = 410.5 cm -3 

The cross section for ([1]) can be determined 
from the Breit-Wigner formula using the measured 
data for the expected resonances. The relativistic 
form can be found in most introductory particle 
physics texts. 



a(s) = a n 
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The peak of the resonance is scaled by <j ma x 
which may be measured experimentally. Mq is the 
mass of the resonance, and T is the width. 

The mean attenuation length for a cosmic 
ray proton of a given energy is the energy fraction 
of scattering times the attenuation length. 
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For photopion production the energy frac- 
tion can be deduced from kinematics, and averages 
roughly AE / E = 15% in the relevant energy regime. 
We can estimate the mean attenuation length of a 
10 20 eV cosmic ray proton, taking an average value 
of a, 500/i6: 
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A cosmic ray traveling a distance, L, would 
have a survival probability less than P SU rvivai ■ 

Psurvivai(L) = 1.0 - exp (-L/L ) 

Therefore, the flux of ultra high energy cos- 
mic rays should be significantly attenuated on cos- 
mological scales. The product of these photohadron 
interactions is a flux of electrons, photons, protons, 
neutrons and neutrinos. However, only the neutri- 
nos travel without attenuation on these cosmological 
distances. These cosmogenic neutrinos could be de- 
tected by Earth based experiments. 

It is a topical goal of the astrophysics com- 
munity to measure and characterize the flux of ul- 
tra high energy neutrinos. Still, no contemporary 
experiment is capable of measuring the cosmogenic 
neutrino flux, so suitable Monte Carlo event genera- 
tors must be available to support ongoing research. 
Here we present the C++ code named gzkfast as 
one possible avenue for simulation of this process. 

3 Program Operation 

GZKFast is principally comprised of two compo- 
nents, a reusable library and a command line 
event generator. The reusable program library, 
libgzkparticle , is written in C++ and designed 
to be modular and portable across platforms. The 
simple invocation program, gzkfast, gives com- 
mand line access to the key functionality provided 
in the GZK library. 

When a user invokes the gzkfast pro- 
gram, a universe is created with a simplified par- 
ticle physics model. Cosmic ray point sources 
are inserted at random locations and managed by 
an EventGenerator thread. The event generator 
thread iterates through the point sources and in- 
jects subsequent ultra high energy particles into 
an EvolutionThread with momentum oriented to- 
wards the Earth. 

Cosmic ray events are managed in parti- 
cle queues of the evolution thread modules using a 
"round robin" strategy. As a result execution time is 
divided relatively equally between available threads. 

Each evolution thread maintains a list of 
relevant Space objects, including the CMB , and ex- 
tragalactic space, or BFieldSpace . The evolution 
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threads also maintain a list of suitable Earth based 
Detector objects. 

Particle evolution continues as particles 
propagate through each of the known spaces and are 
given a chance to be detected by known detectors. 
Basic particle kinematics are used for propagation. 
In the extragalactic B field, assumed to be static and 
uniform over distance dx, the kinematics of particle 
propagation are governed by a few simple relations. 

dp e 

— = -v x B 

dt c 

d Pr 

p = p+-5t 
x = x + vSt 
t = t + St 

The time step, St, is defined based on the 
user specified distance step, dx. Particle interactions 
with CMB photons are evaluated with an accept- 
reject strategy. Then after propagation, each parti- 
cle is allowed to decay with probability Pdecay 

Pdecay(t) = 1.0 - exp (~*/t) 

Decay products are reinserted into the evo- 
lution thread particle queue. If no decay or parti- 
cle production takes place during propagation, each 
known detector is tested to see if it can see the par- 
ticle. If the particle "hits" the detector the particle 
becomes an "event." 

Each event is logged to a relevant output 
location. Particle evolution continues until the user 
specified number of events have been detected. Once 
the specified number of events have been recorded, 
event generation is halted. Program execution then 
terminates once all particle queues are empty. 

4 Basic Results and Spectra 

Each run of gzkfast is configurable to allow the 
user to explore different aspects of cosmogenic neu- 
trino phenomenon. With suitable configuration it is 
easy to produce a wide variety of relevant data. 

GZKFast uses integrated data provided by 
the Particle Data Group [7] to sample the proton- 
photon cross section versus y/s. Accelerator data 
provided in the PDG dataset are linearly interpo- 
lated from point to point to produce a complete dis- 
tribution as depicted in Figure 2] 



Using to determine the cross section for 
the A1232 and Ai6oo- These data are tabulated 
in the file specified by the "cmb" basename, with 
"_sigmaplot.dat" appended. 

Simulated py Cross Section 
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Figure 1: The ^7 cross section measured in accel- 
erator experiments (Black) [7J. The 1232 MeV A 
resonance as determined from the Breit-Wigner for- 
mula (Red). The difference is taken to be the cross 
section for direct pion photoproduction. 

The cosmic microwave background is uni- 
formly sampled as a black body distribution at 
2.725 K [6]. These data are output in the file speci- 
fied by the "cmb" basename, with "_photonhist.dat" 
appended. 

The Number Distribution of the CMB 




0.003 



Energy [eVl 

Figure 2: The number of photons at a given energy 
for a 5000 v event simulation. 

At program termination the gzkfast pro- 
gram produces a variety of output spectra for the 
various processes that it tracked. Most importantly 
GZKFast is capable of producing neutrino spec- 
tra resulting from super-GZK cosmic rays. These 
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data are output in the file specified by the "neu- 
trino" basename, with "_hist.dat" or "_event.dat" 
appended depending on if the data set is a histogram 
or event file. 



Neutrino Flux from CRs 
10 i i i i i 1 1 <m • • 




Figure 3: The neutrino histogram gives the energy 
distribution of an expected neutrino flux. 

GZKFast also provides individual event 
data for protons and neutrons, neutrinos and pho- 
tons. These event data can be utilized in detector 
simulation or since they also include right ascension 
and declination information they can be used to pro- 
duce sky maps with variable sources, distances, and 
other configuration differences. 



GZKFast source code is available to the public upon request. 
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6 GZKFast Reference 
6.1 Class Listing 

CMB The Cosmic Microwave Background Object samples the cross section of pj and 

the distribution of a user specified black body spectrum to determine if an 
interaction with a cosmic ray will occur. 

CMBDist Sample the distribution of a black body radiation of a given temperature. 

Delta Program representation of a A (1232 MeV) particle. Provides decay kinematics 

and particle properties. 

Deltal600 Program representation of a A (1600 MeV) particle. Provides decay kinematics 

and particle properties. 

Detector Abstract virtual class interface for GZKFast spherical detectors. A detector is 

a volume of space which can be hit by a simulated particle. 

Electron Program representation of an e~ particle. Provides decay kinematics and par- 

ticle properties. 

EventGenerator Threaded object which iterates through cosmic ray sources producing parti- 
cle events from a predetermined energy spectrum and inserting them into the 
particle queues. 

EvolutionThread Threaded object which iterates through a specific particle event queue and 
propagates the particles through one of an arbitrary number of spaces. After 
the particle propagates, the detectors are given a chance to detect it. 

G4Vector A 4 vector representation. 

GGuard Guard a critical section. 
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GHistogram 

GLogHistogram 

GMath 

GMatrix 

GMutex 

GRunThread 

GThread 

GVector 

GVegas 

BFieldSpace 

Mu 

Neutron 
Nu 

NuDetector 

Particle 

Photon 

PhotonDetector 
Pion 

PiZero 

Proton 

ProtonDetector 

ProtonSource 

ProtonSpectrum 

Source 

Space 

Sphere 

ThreeBodyDecay 
TwoBodyDecay 

Universe 



A linear scale histogram class. 
A log scale histogram class. 

Provides rudimentary numerical operations and random number generation. 

An object representation of an n x n matrix. 

An object to provide mutual exclusion across platforms. 

A CH — h implementation of a Java-like "java.lang. Thread" object. 

An object to provide rudimentary threading support across platforms. 

A vector of arbitrary size, based on std: :valarray. 

The Vegas Numerical Recipe. 

Propagate a particle through a "static" magnetic field using Hamiltonian for- 
malism. 

Program representation of an /x particle. Provides decay kinematics and particle 
properties. 

Program representation of an neutron. Provides decay kinematics and particle 
properties. 

Program representation of a neutrino. Provides decay kinematics and particle 
properties. 

An instance of Detector which is capable of detecting neutrinos. 
An abstract virtual class interface for particle objects. 

Program representation of a photon. Provides decay kinematics and particle 
properties. 

An instance of Detector which is capable of detecting photons. 
Program representation of a ir. Provides decay kinematics and particle proper- 
ties. 

Program representation of a 7r°. Provides decay kinematics and particle prop- 
erties. 

Program representation of a p + . Provides decay kinematics and particle prop- 
erties. 

An instance of Detector which is capable of detecting protons. 
A point source for cosmic ray protons. 

An object to sample a power law energy spectrum of cosmic ray protons. 
An abstract virtual class interface for cosmic ray sources. 

An abstract virtual class interface for spaces. A space allows a particle to 
propagate. 

A geometric representation of a sphere. Provides ray intersection for detector 
operation. 

A generalized Monte Carlo of the reaction A — > B + C + D. 

A solution for two body decay in the ultra high energy limit. Decay products 

are collinear with parent particles. 

Provides a frame for calculating redshift and volume integrals. 



6.2 Command Line Reference 



Parameter 


Units 


Default 


Description 


-sources 


Number 


3 


The number of cosmic ray sources 
to add to the simulation. Sources 
are added randomly to a shell of 
size specified by -near and -far. 


-events 


Number 


100 


The number of neutrino events to 
simulate before stopping. 
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-threads 


Number 


2 


The number of asynchronous exe- 
cution paths simultaneously pro- 
cessing input events. 


-alpha 


Number 


-2.7 


Simulate a proton input spec- 
trum of the form E a . 


-low 


EeV 


50 


The starting energy for the input 
proton distribution. 


-hi 


EeV 


50000 


The ending energy for the input 
proton distribution. 


-near 


Mpc 


150 


The distance to the nearest cos- 
mic ray source. 


-far 


Mpc 


200 


The distance to the furthest cos- 
mic ray source. 


-dx 


Kpc 


250 


The distance corresponding to 
one iteration or step of the Monte 
Carlo integration. 


-rad 


Kpc 


250 


The radius of the "detector" vol- 
ume. The "detector" is defined 
to be a volume of given radius 
centered about the earth. Any 
particle which will intersect the 
"detector" is considered to be an 
"event." 


-bfield 


Gauss 


le-9 


The maximum magnitude of the 
uniform extragalactic magnetic 
field. 


-quality 


Number 


5e-9 


The precision of the beta decay 
Monte Carlo. 


-proton 


File Name 


"proton" 


The base name of the files for 
writing proton events and his- 
togram output. 


-2nd 


File Name 


"secondaryproton" 


The base name of the files for 
writing secondary proton events 
and histogram output. 


-V 


File Name 


"neutrino" 


The base name of the files for 
writing neutrino events and his- 
togram output. 


-cmb 


File Name 


"cmb" 


The base name of the files for 
writing the cmb photon energy 
histogram and cross section sam- 
ple data. 


-photon 


File Name 


"photon" 


The base name of the files for 
writing photon events and his- 
togram output. 



6.3 Program Output 

$ gzkfast: -[ arguments ] ... 

Simulate a flux of ultra high energy neutrinos from cosmic ray sources, 
-sources # - The number of sources. 
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-events # 


- The number of events to simulate. 


—threads # 


— The number of processor threads . 


—alpha $ 


— Simulate E alpha spectrum. 


—low # 


— Least energy [EeV] . 


-hi ff 


- Highest energy [EeV] . 


-near # 


- Least source distance [Mpc] . 


-tar # 


- Highest source distance [Mpc] . 


ax # 


— distance step [Kpc] . 


-rad # 


- detector radius [Kpc] . 


-bfield # 


- B field strength [Gauss] . 


-quality # 


- The precision of Monte Carlo convergence 


-proton file 


- Name of file for input protons. 


-2nd file 


- Name of file for secondary protons. 


-v file 


- Name of file for neutrinos 


-cmb file 


- Name of file for cmb distributions 


-gamma file 


- Name of file for photons. 



6.4 File Formats 

Energy [eV] E dN/dE [cm" 2 s" 1 sr" 1 ] 

4.235324597282258e+16 6.205560549060830e-05 

4.829213438776459e+16 5.812199626341997e-05 

5.506379003919574e+16 1.066115138632034e-04 

Table 3: Example neutrino_hist.dat output. 



RA [deg] Dec [deg] 

73.8507 -31.1954 

73.8507 -31.1954 

82.8030 -90.6093 

82.8030 -90.6093 



E [eV] p x [eV] 

4.2491+16 2.6999+16 

8.9698+18 5.6996+18 

1.0938+16 -3.6303+15 

1.4671+19 -4.8691+18 



Pv [eV] Pz [eV] 

-7.5372+15 -3.1932+16 

-1.5911+18 -6.7409+18 

3.6691+15 -9.6441+15 

4.9212+18 -1.2934+19 



Table 4: Example neutrino_event.dat output. 



References 



[i] 

[2] 
[3] 
[4] 
[5] 
[6] 
[7] 



A. Miicke, et al.,|astro-ph"/9903478 



K. Greisen, Phys. Rev. Lett. 16, 748 (1966). 

G. T. Zatsepin, V. A. Kuz'min, Zh. Eksperim. i Teor. Fiz. - Pis'ma Redakt. 4, 78 (1966). 



A. Miicke, J.P Rachen, ct al., astro-ph/9808279 



F. W. Stecker, Phys. Rev. Lett. 21, 14 (1968). 
C.L. Bennett et al., ApJS 148, 1 (2003). 

S. Eidelmann et al., Particle Data Group, Phys. Lett. B592, 1 (2004). 



R. Engel, D. Seckel, and T. Stanev, [astro^ph /0101216 



